%
% demo_with_random_vectors: this file demonstrates the 
% usage of the MSCD_Wav toolbox on an example of 
% multi-scale graph. It specifically uses the wavelet 
% transform of several random vectors to compute the 
% partitions and their stability.

% This file is part of the MSCD_Wav toolbox (MultiScale 
% Community Detection using Wavelets toolbox)
% Copyright (C) 2014, Nicolas Tremblay. 
%
% The MSCD_Wav toolbox is free software: you can 
% redistribute it and/or modify it under the terms of 
% the GNU General Public License as published by
% the Free Software Foundation, either version 3 
% of the License, or (at your option) any later version.
%
% The MSCD_Wav toolbox is distributed in the hope that 
% it will be useful, but WITHOUT ANY WARRANTY; without 
% even the implied warranty of MERCHANTABILITY or FITNESS 
% FOR A PARTICULAR PURPOSE.  See the GNU General Public 
% License for more details.
%
% You should have received a copy of the GNU General 
% Public License along with the MSCD_Wav toolbox.  
% If not, see <http://www.gnu.org/licenses/>.

%% load your adjacency matrix (example here: Sales Pardo graph)
clear all; close all;
%S3=10;S2=30;S1=120; N=640;
%COM3=repmat([1:N/S3]',1,S3)'; COM3=COM3(:);
%COM2=repmat([1:N/(S2+S3)]',1,S2+S3)'; COM2=COM2(:);
%COM1=repmat([1:N/(S1+S2+S3)]',1,S1+S2+S3)'; COM1=COM1(:);
%COM=[COM1,COM2,COM3];
%rho=1; kbar=16; % parameters of SP graph
%[A] = SP_hbenchmark(COM, rho, kbar); % creates SP graph
A = load('celegans_matlab_graph.mat').A

A
%% Multi scale analysis (once adjacency matrix A is loaded)

% automatically detect the largest component : %
N=length(A);
[CC] = conn_comp (A,1);
removenodes=[];
if size(CC,1)~=1
    warning(['your graph is disconnected in ' num2str(size(CC,1)) ' components : calculations are done on the largest component. Size of components : ' char(10) '' num2str(sum(CC~=0,2)') ]);
    Aold=A;
    [~,ind]=max(sum(CC~=0,2));
    keepnodes=sort(CC(ind,CC(ind,:)~=0));
    A=A(keepnodes,keepnodes);
    removenodes=setdiff([1:N],keepnodes);
end

Nscales=100; %number of scales
Nsignals=50; % number of random signals used
Jstab=100; %number of random vector sets used to compute stability
lap_flag='normlap';
% equals 'lap' if one wants to use the classical Laplacian L;
% equals 'normlap' if one wants to use the normalized 
% Laplacian D^(-0.5)LD^(-0.5);
% equals 'rwlap' if one wants to use the "random walk" Laplacian D^(-1)L.
filter_flag='scaling_function';
% equals 'wavelet' if one wants to use correlation of wavelets; 
% equals 'scaling_function' if one wants to use correlation of scaling 
% functions. 

display('computing partitions and their stability...');
[tg,COM_RV,stab]=MSCD_partition_RV(Nscales,A,Nsignals,Jstab,lap_flag,filter_flag);
% * tg is the vector of the Nscales scales between 
% the automatically computed scale boundaries
% * COM is a NxNscales matrix where the i-th column 
% corresponds to the community structure found at 
% scale tg(i).
% * stab is the stability measure of each partition.
% stab(i) is the stability of the partition found at scale 
% tg(i).

save('data.mat')

% to visualize partition results in GEPHI (for example), create 
% the network file Partition_RV.gdf
display('creating files containing results...');
MSCD_create_gdf(A,COM_RV,'Partition_RV.gdf');

%% plot instability:
f1=figure; hold on;
plot(log10(tg),1-stab,'linewidth',2);
set(gca,'fontsize',20);
set(gca,'ylim',[0,1]);
set(gca,'xlim',[min(log10(tg)),max(log10(tg))]);
aux=get(gca,'xtick');set(gca,'xticklabel',{floor((10.^aux).*10)/10});
xlabel('scale s');
title('Instability')
ylabel('1-\gamma_a');
set(f1,'position',[444   450   560   322]);
set(gcf,'PaperPositionMode','auto');

%% one can compute and plot performance if ground truth exists:
% (here ground thruth is represented by COM1 COM2 and COM3)
for t=1:Nscales
    isCOM1_RV(t)=PartAgreeCoef_ARonly(COM_RV(:,t),COM1);
    isCOM2_RV(t)=PartAgreeCoef_ARonly(COM_RV(:,t),COM2);
    isCOM3_RV(t)=PartAgreeCoef_ARonly(COM_RV(:,t),COM3);   
end

f1=figure; hold on;
plot(log10(tg),isCOM1_RV,'rs','linewidth',2);
plot(log10(tg),isCOM2_RV,'b<','linewidth',2);
plot(log10(tg),isCOM3_RV,'ko','linewidth',2);
legend('Large Scale','Medium Scale','Small Scale');
set(gca,'fontsize',20);
set(gca,'ylim',[0,1]);
set(gca,'xlim',[min(log10(tg)),max(log10(tg))]);
aux=get(gca,'xtick');set(gca,'xticklabel',{floor((10.^aux).*10)/10});
xlabel('scale s');
title('Multiscale Community Mining Results')
ylabel('Adj. Rand index');
set(f1,'position',[444   450   560   322]);
set(gcf,'PaperPositionMode','auto');
